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Abstract 

Using a simple, annealed model, some of the key features of the recently introduced extremal 
optimization heuristic are demonstrated. In particular, it is shown that the dynamics of local 
search possesses a generic critical point under the variation of its sole parameter, separating phases 
of too greedy (non-ergodic, jammed) and too random (ergodic) exploration. Comparison of various 
local search methods within this model suggests that the existence of the critical point is essential 
for the optimal performance of the heuristic. 
PACS number(s): 02.60.Pn, 05.40.-a, 64.60. Cn, 75.10.Nr. 
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I. INTRODUCTION 



Many situations in physics and beyond require the solution of NP-hard optimization 
problems, for which the typical time needed to ascertain the exact solution apparently grows 
faster than any power of the system size [1] . Examples in the sciences are the determination 
of ground states for disordered magnets [2, 3, 4, 5, 6, 7] or of optimal arrangements of atoms 
in a compound [8] or a polymer [9, 10, 11]. With the advent of ever faster computers, the 
exact study of such problems has become feasible [12, 13]. Yet, with typically exponential 
complexity of these problems, many questions regarding those systems still are only acces- 
sible via approximate, heuristic methods [14] . Heuristics trade off the certainty of an exact 
result against finding optimal or near-optimal solutions with high probability in polynomial 
time. Many of these heuristics have been inspired by physical optimization processes, for 
instance, simulated annealing [15, 16] or genetic algorithms [17]. 

Extremal optimization (EO) was proposed recently [2, 18], and has been used to treat a 
variety of combinatorial [19, 20, 21] and physical optimization problems [7, 22, 23]. Com- 
parative studies with simulated anneahng [18, 19, 24] and other Metropolis [25] based 
heuristics [26, 27] have established EO as a successful alternative for the study of NP- 
hard problems. EO has found a large number of applications, for instance, in pattern 
recognition [28, 29], signal filtering [30], transport problems [31] molecular dynamics sim- 
ulations [32], artificial intelligence [33, 34], social modeling [35], and 3d spin glass mod- 
els [26, 36, 37]. There are also a number of studies that have explored basic features of the 
algorithm [27, 38], extensions [39, 40, 41], and rigorous performance properties [42, 43]. 

In this article, we will use a simple, annealed model of a generic combinatorial opti- 
mization problem, introduced in Ref. [38], to compare analytically certain variations of local 
search with EO, and of Metropolis algorithms such as Simulated Annealing (SA) [15, 16]. 
This comparison affirms the notion of "optimization at the ergodic edge" that motivated the 
T-EO implementation [7, 18]. This implementation possesses a single tunable parameter, 
T, which separates a phase of "greedy" search from a phase of wide fluctuations, combin- 
ing both features at the phase transition into an ideal search heuristic for rugged energy 
landscapes [44]. The model helps to identify the distinct characteristics of different search 
heuristics commonly observed in real optimization problems. In particular, revisiting the 
model with a "jammed" state from Ref. [38] proves the existence of the phase transition to 
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be essential for the superiority of EO, at least within a one-parameter class of local search 
heuristics. At the phase boundary, EO descends sufficiently fast to the ground state with 
enough fluctuations to escape jams. 

This article is organized as follows: In the next section, we introduce the annealed opti- 
mization model, followed in Sec. Ill by a short review of the local search heuristics studied 
here, in particular, of one-parameter variations of EO and of Metropolis-based search. Then, 
in Sec. IV we compare our analytical results for each heuristic in the annealed model. In 
Sec. V we show why versions of EO lacking a phase transition fail to optimize well. We 
summarize our results and draw some conclusions in Sec. VI. 

II. ANNEALED OPTIMIZATION MODEL 

As described in Ref. [38], we can abstract certain combinatorial optimization problems 
into a simple, analytically tractable model. To motivate this model we imagine a generic 
optimization problem as consisting of a number of variables 1 < i < n, each of which 
contributes an amount — Aj to the overall cost per variable (or energy density) of the system, 

-~Ea.. (1) 

i=l 

(The factor 1/2 arises because local cost are typically equally shared between neighboring 
variables.) We call A, < the "fitness" of the variable, where larger values are better and 
Aj = is optimal for each variable. Correspondingly, e = is the (optimal) ground state 
of the system. In a realistic problem, variables are correlated such that not all of them 
could be simultaneously of optimal fitness. But in our annealed model, those correlations 
are neglected. 

A concrete example for the above definitions is provided by a spin glass with the Hamil- 
tonian 

if = - ^ Jijaiaj (2) 

with some quenched random mix of bonds Jjj =e {—1, 0, -1-1} and spin variables = ±1 [7]. 
With Aj = (Tj JijCj — \Jij\, counting (minus) the number of violated bonds of each 
spin i (among its ojj non-zero bonds), it is e = H/n + eo, where eo is an insignificant constant. 
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We will consider that each variable i is in one of a + 1 {ai — a constant) different fitness 
states Xi. We can specify occupation numbers Ua, < a < a, for each state a, and define 
occupation densities pa = ria/n (a = 0, . . . , a). Hence, any local search procedure [14] with 
single- variable updates, say, can be cast simply as a set of evolution equations for the Pa{t), 
i. e. 

Pb^^^Th^aQa- (3) 

a 

The Qa are the probabilities that a variable in state a gets updated; any local search process 
(based on updating a finite number of variables) defines a unique set of Qa-, as we will see 
below. The matrix T^^a specifies the net transition to state b given that a variable in state a is 
updated. This matrix allows us to design arbitrary, albeit annealed, optimization problems. 
Both, T and Q generally depend on the Pa{t) as well as on i explicitly. 

We want to consider the different fitness states equally spaced, as in the spin glass example 
above, where variables in state a contribute aI\E to the energy to the system. Here AE > 
is an arbitrary energy scale. Thus minimizing the "energy" density 

^^^E^-^^^O' (4) 

a 

defines the optimization problem in this model. Conservation of probability and of variables 
implies the constraints 

J2Qa = h J2'^a,b = {0<b<a), (5) 

a a 
a a 

While this annealed model eliminates most of the relevant features of a truly hard opti- 
mization problem, such as quenched randomness and frustration [45], two basic features of 
the evolution equations in Eq. (3) remain appealing: (1) The behavior of a system with a 
large number of variables can be abstracted into a relatively simple set of equations, describ- 
ing their dynamics with a small set of unknowns, and (2) the separation of update process, 
T, and update preference, Q, lends itself to an analytical comparison between different 
heuristics. 
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III. LOCAL SEARCH HEURISTICS 



The annealed optimization model is quite generic for a class of combinatorial optimization 
problems. But it was designed in particular to analyze the "Extremal Optimization" (EO) 
heuristic [38], which we will review next. Then we will present the update probabilities Q 
through which each local search heuristic enters into the annealed model in Sec. II. Finally, 
we also specify the update probabilities Q for Metropolis-based local searches, such as SA. 

A. Extremal Optimization Algorithm 

Here we only give a quick review of the EO heuristic as we will use it below. More 
substantive discussions of EO can be found elsewhere [2, 7, 18]. EO is simply implemented 
as follows: For a given configuration {ai}f^i, assign to each variable (jj an "fitness" 

A, = -0,-1, -2,..., -a (7) 

(e. g. Aj = —{^violated bonds} in the spin glass), so that Eq. (1) is satisfied. Each variable 
falls into one of only a + 1 possible states. Say, currently there are variables with the 
worst fitness, A = —a, Ua-i with A = — (a — 1), and so on up to no variables with the 
best fitness A = 0. (Note that n = ^j^j.) Select an integer k {1 < k < n) from some 
distribution, preferably with a bias towards lower values of k. Determine < a < a such 
that X^6=(j+i JT-b < k < Ylb=a''^b- Notc that lower values of k would select a "pool" Ua with 
larger value of a, containing variables of lower fitness. Finally, select one of the Ua variables 
in state a and update it unconditionally. As a result, it and its neighboring variables change 
their fitness. After all the effected A's and n's are reevaluated, the next variables is chosen 
for an update, and the process is repeated. The process would continue to evolve, unless an 
extraneous stopping condition is imposed, such as a fixed number of updates. The output 
of local search with EO is the best configuration, with the lowest e in Eq. (1), found up to 
the current update step. 

Clearly, a random selection of variables for such an update, without further input of 
information, would not advance the local search towards lower-cost states. Thus, in the 
"basic" version of EO [18], each update one variable among those of worst fitness would be 
made to change state (typically chosen at random, if there is more than one such variable). 
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This provides a parameter-free local search of some capability. But variants of this basic 
ehmination-of-the-worst are easily conceived. In particular, Ref. [18] already proposed r-EO, 
a one-parameter (r) selection with a bias for selecting variables of poor fitness on a slowly 

varying (power-law) scale over the ranking 1 < k < n of the variables by their Aj. In detail, 
r-EO is characterized by a power-law distribution over the fitness-ranks k, 

Pr{k) ^ ~ I ^ k-^ {l<k<n). (8) 
1 — ^ 

It is a major point of this paper to demonstrate the usefulness of this choice. Hence, we 
will compare the effect of this choice with a plausible alternatives, //-EO, which uses an 
exponential scale, 

^.(fc) = (l<fc<n). (9) 

In fact, we show that the exponential cut-off in //-EO, which is fixed during a run, provides 
inferior results to r-EO. Unhke r-EO, ji-EO does not have a critical point affecting the 
behavior of the local search. 

Although Ref. [42] has shown rigorously, that an optimal choice is given by using a sharp 
threshold when selecting ranks, the actual value of this threshold at any point in time is 
typically not obvious (see also Ref. [43]). We will simulate a sharp threshold s {1 < s <n) 
via 

for r — > oo. Since we can only consider fixed thresholds s, which gives results similar in 
character to //-EO, it is not apparent how to shape the rigorous results into a successful 
algorithm. 



B. Update Probabilities for Extremal Optimization 

As described in Sec. Ill A (and in Ref. [38]), each update of r-EO a variable is selected 
based on its rank according to the probability distribution in Eq. (8). When a rank k{< n) 
has been chosen, a variable is randomly picked from state a, if k/n < pa, from state a — 1, 
if Pa < k/n < Pa + Pa-i, ^-^d so on. We introduce a new, continuous variable x — k / n, for 
large n approximate sums by integrals, and rewrite P{k) in Eq. (8) as 

p(^) = (-<x<i), (11) 
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where the maintenance of the low-x cut-off at 1/n will turn out to be crucial. Now, the 
average likelihood in EO that a variable in a given state is updated is given by 

Qa = p{x)dx = — (p^"^ - n^-^) , 

Ji/n l-n^ ^ 

rpa+pa-1 

Qa-1 = j p{x)dx = ^ _ + PoS"'^ - Pa~^] > 

Pol 



Qo = r pix)dx ^ [1 - {1 - po)'-^] , (12) 
Ji-po 1 - n 

where in the last hne the norm Yla Pa — ^ was used. These values of the Q's completely 
describe the update preferences for r-EO at arbitrary r. 

Alternatively, if we consider the //-EO algorithm introduced in Eq. (9), we have to replace 
the power- law distribution in Eq. (11) with an exponential distribution: 

Hence, for p-EO we have 

X _ g-M(Pa-l/n) 
^ 1 _ e-/.(l-l/n) ' 

Q-P'iPa-i/n) _ g-M(pa+Pa_l-l/n) 

^"-1 ^ 1 _ g-Ml-l/n) ' 

g-^(l-po-l/n) _ g-/i(l-l/n) 

^0 = 1 - e-Mi-Vn) • (14) 

Similarly, we can proceed with the threshold distribution in Eq. (10) to obtain 

Psix) OC 7 7 (— < X < 1), (15) 

' ]^ _|_ gr(nx-s) n~ ~ 

with some proper normalization. While all the integrals to obtain Q are elementary, we do 
not display the rather lengthy results here. 

Note that all the update probabilities in each variant of EO are independent of T (i. e. 
any particular model), which remain to be specified. This is quite special, as the following 
case of Metropolis algorithms shows. 



C. Update Probabilities for Metropolis Algorithms 



It is more difficult to construct Q for Metropolis-based algorithms [25] hke simulated 
annealing [15, 16]. Let's assume that we consider a variable in state a for an update. 
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Certainly, Qa would be proportional to pa, since variables are randomly selected for an 
update. The Boltzmann factor e'^^^" for the potential update from time i — > i + 1 of a 
variable in a, aside from the inverse temperature Pit), only depends on the entries for T^^: 



n 
2 



n 



n 
2 
n 
2 



Y^^Pbit) 
b 

b c 



(16) 



where the subscript a expresses the fact that it is a given that a variable in state a is 
considered for an update. Hence, we find for the average probability of an update of a 
variable in state a 



Qa = -^PaHiin <( l,exp 



(17) 



where the norm M is determined via YlaQa ~ ^- Unlike for EO, the update probabilities 
for SA are model-specific, i. e. depend on T. 



IV. COMPARISON OF LOCAL SEARCH HEURISTICS 

To demonstrate the use of these equations, we consider a simple model of an energetic bar- 
rier with only three states {a — 2) and a constant flow matrix Th^a — [—Sb,a + S{2+b mod 3),a\/'n, 
depicted in Fig. 1. Here, variables in pi can only reach their lowest-energy state in po by 
first jumping up in energy to p2- Eq- (3) gives 

Po = - (-Qo + Q2) , Pi = -{Qq-Qi), P2 = -(Qi-Q2), (18) 
n n n 

with Q discussed in Sec. IIIB for the variants of EO. 

Given T, we can now also determine the update probabilities for Metropolis according 
to Eqs. (17). Note that for a = 2 we can evaluate the min as 1, since "^j^bTh^a < always. 
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Flow up 

FIG. 1: Flow diagram with energetic barri- 
ers. Arrows indicate the average number of 
variables transfered, nTh^a, from a state a to 
a state b, given that a variable in a gets up- 
dated. Diagonal elements Ta,a correspond- 
ingly are negative, accounting for the out- 
flow. Note that variables transferring from 
pi to po most first jump up in energy to p2- 




Flow jam 

FIG. 2: Same as Fig. 1, but with a model 
leading to a jam. Variables can only transfer 
from p2 to Po through pi, but only if pi < 9. 
Once pi = 9, flow down from p2 ceases until 
pi reduces again. 



while for a = 0,1 the min always evaluates to the exponential. Properly normalized, we 
obtain 



Qo 



-/3/2 



(l-e-/^/2)p2 + e-^/2' 

P2 



pie 



-0/2 



(1 - e-/3/2)p2 + e-/3/2 



(19) 



(l_e-/3/2)p2 + e-^/2' 

It is now very simple to obtain the stationary solution: For p = 0, Eqs. (18) yields 
Qf) = Qi = Q2 = l/^i and we obtain from Eq. (12) for r-EO: 



Po = 1 
for p-EO: 



1 r-l — 

-rf ^ + -\ , p2 



2 1^ 



pi = 1 - Po - P2, 



(20) 



Po = - In 
P 



3^3 



P2 = 1 In 

P 



1 + 2 

3^3 



pi = 1 - Po - P2, (21) 



and for Metropolis: 



For EO with threshold updating, we obtain 
^° 3 3n n 3nr ^ 



1 



(e"" + e''") (1 + e"(^-')) ^ + e^^^^+i) (i + e'-("-^)) ^ 



H In 

nr 

3 3n n 3nr 



H In 



and, assuming a threshold anywhere between l<s<n, forr— >oo: 

^ 2s + 1 s + 2 ^ 

Po = 1 P2 = Pi = 1 - Po - P2 (24) 

Therefore, according to Eq. (4), Metropolis reaches its best, albeit sub-optimal, cost 
e = l/4>0at/9^oo, due to the energetic barrier faced by the variables in pi, see Fig. 1. 
(Since fluctuations from the mean are suppressed in this model, even a slowly decreasing 
temperature schedule as in Simulated Annealing would not improve results.) In turn, //-EO 
does reach optimality (po = 1, hence e = 0), but only for p — > oo. Note that in this limit, 
p-EO reduces back to the "basic" version of EO discussed in Sec. Ill A. The result for 
threshold updating in EO are more promising: near-optimal results are obtained, to within 
0{l/n), for any finite threshold s. But again, results are best for small s — > 1, in which 
limit we revert back to "basic" EO. 

The result for r-EO is most remarkable: For n — > oo at r < 1 EO remains sub-optimal, 
but reaches the optimal cost for allr > 1! As discussed in Ref. [38], this transition at r = 1 
separates an (ergodic) random walk phase with too much fluctuation, and a greedy descent 
phase with too little fluctuation, which would trap r-EO in problems with broken ergodicity 
[46]. This transition derives generically from the scale-free power-law in Eq. (8), as was 
already argued on the basis of numerical results for real NP-hard problems in Refs. [18, 20]. 

V. JAMMING MODEL FOR /x-EO 

In this section, we revisit the "jammed" model treated in Ref. [38] for r-EO and repeat 
that calculation for p-EO. As in the example in Sec. IV, p-EO proves inferior to r-EO: 
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Lacking the phase of optimal performance in the r-parameter space, the required fine-tuning 
of II does not succeed in satisfying the conflicting constraints imposed on the search. 

Naturally, the range of phenomena found in a local search of NP-hard problems is not 
limited to energetic barriers. After all, so far we have only considered constant entries for 
Th^a- Therefore, in our next model we want to consider the case of T depending linearly on 
the Pi discussed in Ref. [38] for r-EO. This model highlights significant differences between 
the T-EO and the /i-EO implementation. 

Prom Pig. 2, we can read off T and obtain for Eq. (3): 



Po 
Pi 



1 
n 
1 
n 



-Qo + 2<5i 



:Qo -Qi + {e- pi)Q2 



(25) 



. and p2 = — Po — Pi from Eq.(6). Aside from the dependence of T on pi, we have also 
introduced the threshold parameter 6. In fact, if ^ > 1, the model behaves effectively like 
the previous model, and for ^ < there can be no flow from state 2 to the lower states 
at all. The interesting regime is the case < ^ < 1, where further flow from state 2 into 
state 1 can be blocked for increasing pi, providing a negative feed-back to the system. In 
effect, the model is capable of exhibiting a "jam" as observed in many models of glassy 
dynamics [47, 48, 49], and which is certainly an aspect of local search processes. Indeed, the 
emergence of such a "jam" is characteristic of the low-temperature properties of spin glasses 
and real optimization problems: After many update steps most variables freeze into a near- 
perfect local arrangement and resist further change, while a finite fraction remains frustrated 
(temporarily in this model, permanently in real problems) in a poor local arrangement [50] . 
More and more of the frozen variables have to be dislodged collectively to accommodate 
the frustrated variables before the system as a whole can improve its state. In this highly 
correlated state, frozen variables block the progression of frustrated variables, and a jam 
emerges. 

Inserting the set of Eqs. (14) for o; = 2 into the model in Eqs. (25), we obtain 



.1 1 



Q 1 
I _ _g/iPO _| gM(pO+Pl) 

2 2 



1 



Pi 



1 3 

I gMPO _ gM(PO+Pl) -\- (^0 — pj^j (g/^Cl-l/^^) _ gMPO+pi)^ 



, (26) 



n e''(i-V'») - 1 

At large times t, the steady state solution, p — 0, yields for po after eliminating pi the 
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implicit equation 

(eM(i-i/n) _ 3gMPo^ ^ (27) 



= ?-?e- + 



^- - ln(3-2e-''''°) 



and according to Eq. (4), again eliminating pi and p2 in favor of po, we can express the cost 
per variable as 

e = l-po-^ln(3-2e-'^''°), 



~ — In 



{p » 1), (28) 



Unlike the corresponding equations in Ref . [38] , which had a phase transition similar to the 
solution for r-EO in Sec. IV, Eqs. (27-28) have no distinct features. In fact, as shown in 
Fig. 3, e{p) behaves similar to the solution for p-EO in Sec. IV: The relation is independent 
of n to leading order and only for ^ ^ oo, po ~^ 1 ^-nd e ^ 0. 

While the steady state (t oo) features of this model do not seem to be much different 
from the model in Sec. IV, the dynamics at intermediate times t is more subtle. In particular, 
as was shown in Ref. [38], a "jam" in the flow of variables towards better fitness may ensue 
under certain circumstances. The emergence of the jam depends on initial conditions, and 
its duration will prove to get longer for larger values of p. If the initial conditions place 
a fraction po > 1 — 9 already into the lowest state, most likely no jam will emerge, since 
Pi{t) < 6 for all times, and the ground state is reached in t = 0{n) steps. But if initially 
Pi + P2 = 1 — Po > ^) and p is sufficiently large, p-EO will drive the system to a situation 
where pi ~ ^ by preferentially transferring variables from p2 to pi. Then, further evolution 
becomes extremely slow, delayed by the //-dependent, small probability that a variable in 
state 1 is updated ahead of all variables in state 2. 

Clearly, this jam is not a steady state solution of Eq. (26). It is not even a meta-stable 
solution since there are no energetic barriers. For instance, simulated annealing at zero 
temperature would easily find the solution m. t — 0{n) without experiencing a jam. In 
reality, a hard problem would most certainly contain combinations of jams, barriers, and 
possibly other features. 

To analyze the jam, we consider initial conditions leading to a jam, pi(0) + P2(0) > 6 and 
make the Ansatz 



p,{t)^e-r){t) 
12 



(29) 




FIG. 3: Plot of the energy (e) averaged over many ^-EO runs after tmax = lOOn updates with 
different initial conditions as a function of for n = 10, 100, 1000, and 10000 and 9 = 1/2. Plotted 
also is e = Yla ^Paf^ ^ a function of fi resulting from the jam-free steady-state solution {t — oo) of 
Eqs. (27-28) n = oo. The plots show little variability with system size for large n, and remain quite 
sub-optimal for finite fi. As for r ^ oo in Ref. [38], (e) ^ 2/3 - 6/2 - 9^/6 = 19/48 fii 0.396 [51] 
for ij, oo. 

with ?7 <^ 1 for i ^ tjam, where tjam is the time at which po 1- To determine tjam, wc apply 
Eq. (29) to the evolution equations in (26) to get 



1 1 

Po 



n e 



H{l-l/n) _ I 



Q 1 

2 2 



(30) 



where the relation for pi merely yields a self-consistent equation to determine sub-leading 
corrections. 

We can now integrate Eq. (30) from t = (assuming that any jam emerges almost 
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instantly) up to ijam, where po — 1: 

The integral is easily evaluated, and we find for large values of /i: 

tjam ^e'^(^-'^-''°(°)) (/.»!). (32) 

Instead of repeating the lengthy calculation in Ref. [38] for the ground state energy averaged 
over all possible initial conditions for finite runtime tmax oc n, we can content ourselves here 
with a few obvious remarks: A finite fraction of the initial conditions will lead to a jam, hence 
will require a runtime tmax ^ ^jam to reach optimality. Yet, to reach a quality minimum, 
say, e ~ 1/n, would require // ~ n ^ 1 according to Eq. (28). Thus, the require runtime to 
resolve the jam would grow exponentially with system size n, since from Eq. (32) tjam ~ e"^" 
with c = 1 — 9 — po{0) > 0, by definition of the jam above. 

In conclusion, /i-EO can never quite resolve the confiicting demands of pursuing quality 
ground states with a strong bias for selecting variables of low fitness (i. e. /i ^ 1) and the 
ensuing lack of fluctuations required to break out of a jam, which drives up fjam- Simulations 
of this model with p-EO in Fig. 3 indeed show that the best results for (e) are obtained 
at intermediate values of p, which converge to a large, constant error for increasing n. In 
contrast, r-EO provides a range near Topt — 1 ~ 1/ ln(n) [38] with small enough r to fluctuate 
out of any jam in a time near-linear in n while still attaining optimal results as it does for 
any r > 1, see e. g. Sec. IV. 



VI. CONCLUSION 



We have presented a simple model to analyze the properties of local search heuristics. 
The model with a simple energetic barrier demonstrates the characteristics of a number of 
these hcmistics, whether athermal (EO and its variants) or thermal (Metropolis) [27]. In 
particular, it plausibly describes a number of real phenomena previously observed for r-EO 
in a tractable way. Finally, in a more substantive comparison on a model with jamming, 
the exponential distribution over fltnesses, p-EO proves unable to overcome the conflicting 
constraints of resolving the jam while finding good solutions. This is in stark contrast with 
the identical calculation in Ref. [38] using a scale- free approach with a power- law distribution 
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over fitnesses in r-EO. In this approach, a sharp phase transition emerges generically between 
an expansive but unrefined exploration on one side ("ergodic" phase), and a greedy but 
easily trapped search on the other ( "non-ergodic" phase), with optimal performance near 
the transition. 
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